#------------------------------------------------------------------------------
# Import libraries
#------------------------------------------------------------------------------

rm(list=ls())
library(LalRUtils)
libreq(tidyverse, data.table, zoo, tictoc, 
       fst, fixest, PanelMatch, patchwork,
       rio, magrittr, janitor, did, panelView, 
       ggiplot, tictoc, binsreg, interflex)
set.seed(42)
theme_set(lal_plot_theme())

#------------------------------------------------------------------------------




#------------------------------------------------------------------------------
# Define Paths
#------------------------------------------------------------------------------

# R studio
setwd( dirname(rstudioapi::getActiveDocumentContext()$path) )
# R default : unccoment if you use default R
# setwd(getSrcDirectory(function(){})[1])

#------------------------------------------------------------------------------




#------------------------------------------------------------------------------
# Load Data
#------------------------------------------------------------------------------

vcf <- fread("vcf_data_complete.csv", sep = ",")
setnames( vcf, "d", "D")
vcf_data <- copy( vcf )
gfc <- fread("gfc_dta.csv" )


#------------------------------------------------------------------------------



#------------------------------------------------------------------------------
# Plot
#------------------------------------------------------------------------------
# %%
vcf = vcf_data[year>= 1990]
fifth_sched_pre_2k = c("Andhra Pradesh", "Chhattisgarh", 
                       "Gujarat", "Himachal Pradesh",
                       "Orissa", "Rajasthan", "Madhya Pradesh")
pri = vcf[year<=1999 & state %in% fifth_sched_pre_2k]
pri[, first_panch_elec := fcase(
  state == "Andhra Pradesh", 1995,
  state == "Chhattisgarh", 1995,
  state == "Gujarat", 1995,
  state == "Madhya Pradesh", 1994,
  state == "Himachal Pradesh", 1995,
  state == "Orissa", 1997,
  state == "Rajasthan", 1995
)]

# %% flip treatment status since PRI treatment is PESA control and vice vers
pri[, time_pri := year - first_panch_elec]
pri[, D := (1 - sch) * (year >= first_panch_elec)]
pri[, nsch := (1 - sch)]


# %%
pri[, never_treated := max(D) == 0, cellid]
pri[, pref_bin := ntile(cover_1990, 10)]
fitter = function(cut) {
  m = feols(forest_index ~ D | cellid + styear, 
            data = pri[pref_bin >= cut], cluster = ~blk)
}

# Get models
models_pri <- lapply(1:10, fitter )
# Function 
sum_stats_pri = function( cutoff ){
  
  # Filter data
  dat = pri[ pref_bin >= cutoff ]
  dat[, never_treated := max(D) == 0, cellid]
  
  # Get control and treatment value for pre Y
  ctrl <- dat[ never_treated == 1 & year < first_pesa_exposure, 
               mean( forest_index ) ]
  treat <- dat[ never_treated == 0 & year < first_pesa_exposure, 
                mean( forest_index ) ]
  
  year_min <- min(dat$year)
  year_max <- max(dat$year)
  # return vector of values
  return( c( ctrl, treat, year_min, year_max ) )
}

# Generate sum statistics for vcfevsamp
controls_pri <- c()
treated_pri <- c()
min_years_pri <- c()
max_years_pri <- c()
vec_cutoff <- c(1:10)
for ( i in 1:length( vec_cutoff ) ){
  res <- sum_stats_pri( vec_cutoff[ i ] )
  controls_pri[ i ] <- res[ 1 ]
  treated_pri[ i ] <- res[ 2 ]
  min_years_pri[ i ] <- res[ 3 ]
  max_years_pri[ i ] <- res[ 4 ]
}

# All models
treatmap =c(
  "D_fra" = "PESA $\\times$ Scheduled $\\times$ FRA ",
  # GFC
  "def_ha"       = "Annual Deforestation in Hectares",
  "village"      = "Village",
  "village[t]"   = "Village + Village TT",
  # VCF
  "forest_index" = "Forest cover index",
  "green_index"  = "Non-forest green index",
  "built_index"  = "Non-forest index",
  "cellid"       = "Pixel",
  "cellid[t]"    = "pixel + pixel TT",
  # both
  "yr"           = "Year",
  "year"         = "Year",
  "styear"       = "State $\\times$ Year",
  "D"            = "PESA $\\times$ Scheduled",
  "block"        = "Block",
  "blk"          = "Block"
)
fitstat_register("n_new", function(x) summary( x )$nobs , 
                 "\\# Observations")
desc = "Deforestation and Forest cover index"
fn = "figure5_regs"; lab = "tab:figure5_regs";
etable(models_pri,
       style.tex = style.tex(main = "base", 
                             depvar.title = "", 
                             model.title = "", 
                             var.title = "\\midrule", 
                             slopes.title = "\\midrule \\emph{Time Trends}", 
                             yesNo = c("$\\checkmark$", "")),
       signif.code = NA,
       fixef_sizes = T, 
       fixef_sizes.simplify = F,
       fitstat = c( "n_new", "r2" ), 
       tex = TRUE,
       dict = treatmap,
       label = lab,
       title = glue::glue("Effects of Forest Rights Act on Forest Index (pri data)"),
       extralines=list(
         "\\midrule \\emph{Summary Statistics}" = c( rep( "", 10 ) ),
         "Mean Pre-Y (Non-Sch)"= c( controls_pri ),
         "Mean Pre-Y (Sch )"   = c( treated_pri ),
         "Dataset"     = c( rep( "VCF", 10 ) ),
         "Timespan"    = c( rep( "1995-2017", 10 ) )
       ),
       file = "main_figure7_tablePanelA.tex", 
       replace = TRUE
)





# %% one-shot treatment in 2008
vcf[, pref_bin := ntile(cover_1990, 10)]
vcf[, time_fra := year - 2008]
vcf[, D_fra := sch * (year >= 2008)]
fitter = function(cut) {
  m = feols(forest_index ~ D_fra | cellid + styear, 
            data = vcf[pref_bin >= cut], 
            cluster = ~blk)
}

# Get models
models_vcf <- lapply(1:10, fitter )

# Function 
sum_stats_vcf = function( cutoff ){
  
  # Filter data
  dat = vcf[ pref_bin >= cutoff ]
  dat[, never_treated := max(D) == 0, cellid]
  
  # Get control and treatment value for pre Y
  ctrl <- dat[ never_treated == 1 & year < first_pesa_exposure, 
               mean( forest_index ) ]
  treat <- dat[ never_treated == 0 & year < first_pesa_exposure, 
                mean( forest_index ) ]
  
  year_min <- min(dat$year)
  year_max <- max(dat$year)
  # return vector of values
  return( c( ctrl, treat, year_min, year_max ) )
}

# Generate sum statistics for vcfevsamp
controls_vcf <- c()
treated_vcf <- c()
min_years_vcf <- c()
max_years_vcf <- c()
vec_cutoff <- c(1:10)
for ( i in 1:length( vec_cutoff ) ){
  res <- sum_stats_vcf( vec_cutoff[ i ] )
  controls_vcf[ i ] <- res[ 1 ]
  treated_vcf[ i ] <- res[ 2 ]
  min_years_vcf[ i ] <- res[ 3 ]
  max_years_vcf[ i ] <- res[ 4 ]
}


# All models
treatmap =c(
  "D_fra" = "PESA $\\times$ Scheduled $\\times$ FRA ",
  # GFC
  "def_ha"       = "Annual Deforestation in Hectares",
  "village"      = "Village",
  "village[t]"   = "Village + Village TT",
  # VCF
  "forest_index" = "Forest cover index",
  "green_index"  = "Non-forest green index",
  "built_index"  = "Non-forest index",
  "cellid"       = "Pixel",
  "cellid[t]"    = "pixel + pixel TT",
  # both
  "yr"           = "Year",
  "year"         = "Year",
  "styear"       = "State $\\times$ Year",
  "D"            = "PESA $\\times$ Scheduled",
  "block"        = "Block",
  "blk"          = "Block"
)
fitstat_register("n_new", function(x) summary( x )$nobs , 
                 "\\# Observations")
desc = "Deforestation and Forest cover index"
fn = "figure7_regs"; lab = "tab:figure7_regs";
etable(models_vcf,
       style.tex = style.tex(main = "base", 
                             depvar.title = "", 
                             model.title = "", 
                             var.title = "\\midrule", 
                             slopes.title = "\\midrule \\emph{Time Trends}", 
                             yesNo = c("$\\checkmark$", "")),
       signif.code = NA,
       fixef_sizes = T, 
       fixef_sizes.simplify = F,
       fitstat = c( "n_new", "r2" ), 
       tex = TRUE,
       dict = treatmap,
       label = lab,
       title = glue::glue("Effects of Forest Rights Act on Forest Index (VCF data)"),
       extralines=list(
         "\\midrule \\emph{Summary Statistics}" = c( rep( "", 10 ) ),
         "Mean Pre-Y (Non-Sch)"= c( controls_vcf ),
         "Mean Pre-Y (Sch )"   = c( treated_vcf ),
         "Dataset"     = c( rep( "VCF", 10 ) ),
         "Timespan"    = c( rep( "1995-2017", 10 ) )
       ),
       file = "main_figure7_tablePanelB.tex", 
       replace = TRUE
)



